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Abstract 

In this work we make a numerical study of the dynamic universality class 
of the Niedermayer algorithm applied to the two-dimensional Potts model with 
2, 3, and 4 states. This algorithm updates clusters of spins and has a free 
parameter, E , which controls the size of these clusters, such that Eq = 1 is 
the Metropolis algorithm and E = regains the Wolff algorithm, for the Potts 
model. For — 1 < Eq < 0, only clusters of equal spins can be formed: we show 
that the mean size of the clusters of (possibly) turned spins initially grows with 
the linear size of the lattice, L, but eventually saturates at a given lattice size 
L, which depends on E . For L > L, the Niedermayer algorithm is in the 
same dynamic universality class of the Metropolis one, i.e, they have the same 
dynamic exponent. For Eq > 0, spins in different states may be added to the 
cluster but the dynamic behavior is less efficient than for the Wolff algorithm 
(Eq = 0). Therefore, our results show that the Wolff algorithm is the best choice 
for Potts models, when compared to the Niedermayer's generalization. 



1. Introduction 

Numerical simulations of physical systems have been used for decades (for 
a review of numerical methods in statistical physics, see Refs. Q and 0). In 
recent years, however, this method has received a renewed interest, due to the 
increase in computational power and, more important, due to the introduction 
of new algorithms. The developments in the algorithms have the goal to allow 
for more efficient simulations, in many different directions: the calculation of the 
density of states using flat histograms, which allows for obtaining information 
at any temperature from one single simulation Q ; the use of bitwise operations 
and storage, which increases by a great amount the speed of the simulation and 
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saves memory [J]; and the introduction of cluster algorithms, which updates 
collections of spins, decreasing the autocorrelation time and reducing critical 
slowing down 0, 0, H, @| ■ 

The Metropolis algorithm Q, for example, which had been the main choice 
of algorithm for a long time, suffers from a severe critical slowing-down effect 
near critical points. This is in part due to the fact that it updates one spin each 
time (therefore, being in the general category of single-spin algorithms). Critical 
slowing down is measured through the dynamic critical exponent z, defined from 
the dependence of the autocorrelation time t on the linear size of the simulated 
lattice, L, at the critical temperature T c . This dependence is assumed to be in 
the form r ~ L z . Therefore, smaller values of z lead to smaller autocorrelation 
times and more efficient algorithms, since more configurations can be used to 
calculate the necessary averages. The Metropolis algorithm, for example, when 
applied to the Ising model in two dimensions, presents z ~ 2.16 

One possible way to overcome the difficulty of critical slowing down is to de- 
sign algorithms which update clusters of spins (the so-called cluster algorithms), 
that may have a much lower value of z: this is the case for the Swendsen-Wang 
[5] and Wolff H a lgorithms, for which z may be zero for the two-dimensional 
Ising model [9l. Il0|. See also Ref. for another possible dependence of the 
autocorrelation time on the lattice size L. An alternative (and generalization) 
to these last two cluster algorithms, the Niedermayer algorithm, was introduced 
some time ago [l2[ but, to the best of our knowledge, has only had his dynamic 
behavior studied in detail for the Ising and XY models [ll|. Our goal in this 
work is to study in a systematic way the dynamic behavior of the Niedermayer 
algorithm, applied to Potts models in two dimensions, for number of states 
q = 2, 3, and 4, for some values of the free parameter Eo (see below), in order to 
determine the best choice of this parameter for these models. Note also that the 



critical temperature for two-dimensional Potts models are exactly known 13 1 
which allows for a more precise determination of z 



The Potts model is defined by the Hamiltonian [13 1 



H = -J § Si,Si, (1) 
<i,j> 

where Si = 1, 2, . . . , q, Vi, the sum is over nearest neighbors on a lattice (in our 
case, a square lattice) and Ss i ,s :j is the Kronecker delta (<5s ij( g, = 1, if Si = Sj, 
and Ssi,Sj = 0, if Si ^ Sj). We treat the ferromagnetic case in this work, i.e, 
J > 0. In two dimensions, the phase transition for this model is a continuous 
one for q < 4. In our study, we restrict ourselves to this interval. 

This work is organized as follows: in the next section we present the Nieder- 
mayer algorithm for the Potts model and its relation to Metropolis' and Wolff's. 
In Section [3] we review some features connected to the autocorrelation time and 
the dynamic exponent z, in Section|4]we present and discuss our results, and in 
the last section we summarize the results. 
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2. The Niedermayer algorithm 



The Niedermayer algorithm was introduced some time ago as a cluster algo- 
rithm, motivated by the possibility to diminish the value of the dynamic critical 
exponent. It may be seen as a generalization of the Wolff or the Swendsen- 
Wang algorithms. The main idea is to build a cluster of spins and accept its 
updating as a single entity, with a parameter E which controls the nature and 
size of the cluster and its acceptance ratio, as explained below. In this work, 
we have chosen to build the clusters according to the Wolff algorithm (they can 
be constructed according to the Swendsen-Wang rule but the results will not 
differ qualitatively in two dimensions and in higher dimensions Wolff algorithm 
is superior to Swendsen- Wang's). 

In the Niedermayer algorithm, a spin in the lattice, which we will call the 
seed spin, is randomly chosen, being the first spin of the cluster. First-neighbors 
of this spin may be considered part of the cluster, with a probability 

P d AE )-{ 1-^-^ , if E ij<Eo 

r add (^j) - j o f otherwise . W 

where K = J/kT, T is the temperature, J is the exchange constant, and £!y 
is the energy between nearest-neighbor spins in unities of J (i.e, Eij = —S SiSj ). 
First-neighbors of added spins may be added to the cluster, according to the 
probability given above. Each spin has more than one chance to be part of the 
cluster, since it may have more than one first-neighbor in it. When no more 
spins can be added, all spins in the cluster have their state changed to a new 
state with an acceptance ratio A. Assuming that, at the frontier of the cluster 
there are m bonds linking spins in the same state in the old configuration and n 
bonds linking spins in the same state in the new configuration (there are also p 
bonds in the border of the cluster which are in different states in the old and in 
the new configurations, but they cancel out in the expression below), A satisfies: 



A(a -> b) 
A(b -> a) 



e K 



1 - Padd(-J) 
1 - Padd(0) 



-\ n—m 



(3) 



where a — >• b represents the possible updating process, from the "old" (a) to 
the "new" (6) state, which differ from the flipping of all spins in the cluster, 
and b — > a represents the opposite move. This expression ensures that detailed 
balance is satisfied 0. 

One has to consider three different intervals for Eq- 

(i) for — 1 < Eo < 0, only spins in the same state as the seed may be added 
to the cluster, with probability P a dd = 1 — e ~~ K ( 1+E °\ The new state of 
the cluster is randomly chosen between the remaining q — 1 states. The 
acceptance ratio (Eq. |3]) cannot be chosen to be one always and is given 
by A = e - KE o( , m-n) ^ if n < m (i_ 6j jf the energy increases when the spins 
in the cluster are changed), or by A = 1, if n > m (i.e, if the energy 
decreases when the spins in the cluster are changed). If Eq = —I, we 



3 



obtain the Metropolis algorithm, since only one-spin clusters are allowed 
(according to Eq. ([2]), P a dd — for E = — 1) and the acceptance ratio is 
A — e ~ KAE for positive AE and 1 otherwise, where AE = (m — n) is the 
difference in energy when the spin is changed, in units of J; 

(ii) for Eq = 0, again only spins in the same state can take part of the cluster, 
with probability P a dd = l — e~ K . Now, the acceptance ratio can be chosen 
to be 1, i.e, the cluster of like spins is always changed. Again, the new 
state of the cluster is randomly chosen between the remaining q — 1 states. 
This is the celebrated Wolff algorithm; 

(iii) for Eq > 0, spins in different states may be part of the cluster. Consider 
a spin already in the cluster: the probability of adding one of its first- 
neighbors to the cluster is P a dd = 1 — e- K ( l + E a) jf they are in the same 
state or P a dd = 1 — e~ KE ° otherwise. The acceptance ratio is again always 
1. To change the state, for each cluster, we randomly choose a Aq between 
f and q — 1 and perform a cyclic sum. Note that for Eq 3> nearly all 
spins will be in the cluster and the algorithm will be clearly inefficient (in 
fact, it will not be ergodic for Eq — > oo). Therefore, we expect that, if the 
optimal choice of Eq is greater than 0, it will not be much greater than 
this value. 

After constructing a cluster, possible updating it and calculating the rele- 
vant thermodynamic functions for the new configuration, the whole process is 
repeated with a new seed spin. In this way, the Markov chain of configurations 
is generated. 

3. Autocorrelation time and dynamic exponent 

To calculate the relevant averages from a numerical simulation, one has to 
build a Markov chain of spin configurations and use data from uncorrelated 
configurations along this chain. Therefore, one important quantity is the auto- 
correlation time for a given quantity, say $(i), obtained from the autocorrelation 
function p(t): 

P (t) = Jim 1 )- <#>][*(t / +t)-<*>]d* / 

= J [$(*')$(*' + <)-<$ > 2 ] dt', (4) 

Since time is a discrete quantity on Monte Carlo simulations, one has to dis- 
cretize the above equation Q: 

P (t) = -J— . J2 mtm'+t)}- 



itmax t) 2 
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It is usually assumed that the autocorrelation function behaves, in its sim- 
plest form, as [2|: 

p(t)=Ae- t /-, (6) 

This hypothesis has to be corroborated by data and, in some cases, more than 
one exponential term is required [l4j . One point worth mentioning is that 
the autocorrelation function is not well behaved for long times, due to bad 
statistics (this is evident from Eq. [5j since few "measurements" are available 
for long times) . Then, one has to choose the region where the straight line will 
be adjusted very carefully and it turns out that the value of r so obtained is 
strongly dependent on this choice. One other possible way to measure r is to 
integrate p(t), assuming a single exponential dependence on (past and forward) 
time, and obtain: 

r = — 

with: 

P(t) 

Eq. when discretized, leads to [lij 

1 

T = — 

2 

The sum in the previous equation cannot be carried out for large values of 
t, since bad statistics would lead to unreliable results. In order to truncate the 
sum at some point, we use a cutoff (see Ref.[15| and references therein), defined 
as the value in time where the noise in the data is clearly greater than the 
signal itself. We then obtain a first estimate of t using Eq. |H] and then make 
the integral of p(t)/p(0) from the value of the cutoff to infinity. A criterion to 
accept the cutoff is that the value of this integral is smaller than the statistical 
uncertainty in calculating r. 

A note on the definition of "time" is worth stressing here. In the Metropolis 
algorithm, time is measured in Monte Carlo steps (MCS) ; one MCS is defined 
as the attempt to flip N spins, where N is the number of spins in the (finite) 
lattice being simulated (in our case, N = L 2 , where L is the linear size of the 
lattice). For cluster algorithms, one unity of time is defined as the "time" taken 
to build and possibly change a cluster, t. A rescaling of the quantity is necessary, 
to be able to compare the results for different values of Eq , namely: 

t MC S=t—jj-, (10) 

where ijwcs is the time measured in MCS and < n > is the mean cluster size. 
Note that, for Metropolis, < n >— 1 and 1 MCS is the "time" taken to try to 
flip N spins, as usual. In this work, this rescaling has been done and all times 
are expressed in MCS. 

Our first attempt was to fit the autocorrelation time to the expected behav- 
ior, namely r ~ L z , in the critical region, where z is the dynamic exponent. We 



3 4^,, 

cP(0) 
,-\t\/r 



(7) 
(8) 
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have taken as the chosen quantity to calculate r the one which better adjusts 
to Eq. [6] Usually, different quantities lead to autocorrelation functions which 
behave differently as a function of time. A typical example is shown in Fig. 
[TJ where both the magnetization and the energy autocorrelation functions for 
the q = 3 Potts model are depicted as functions of time, for the Niedermayer 
algorithm with Eq = —0.25 and linear sizes L = 16 (main graph) and L = 128 
(inset). The magnetization autocorrelation function presents an abrupt drop 
for small times and L = 16. This shows that this function is not properly de- 
scribed by a single exponential. On the other hand, the energy autocorrelation 
time follows a straight line even for the smallest times. Therefore, we should 
calculate r from the latter, for L = 16, using Eq. [5J Note, however, that, when 
L is increased, the picture changes and now the magnetization autocorrelation 
function is well described by a single exponential (for small and intermediate 
values of time), as depicted in the inset of Fig. Q] Whenever a crossover like 
this is present, we measure the dynamic exponent from the behavior for large 
values of L and for the function which is well described by a single exponential 
for this range of L, using Eq. [5J The fact that, for intermediate values of t, the 
slopes of both curves in Fig. Q] (main graph and inset) appear to be the same, 
is an indication that the autocorrelation times for both are the same. However, 
we have already commented on the drawback of calculating r from the slope of 
the autocorrelation function on a semi-log graph. As final note, we would like 
to mention that we used helical boundary conditions and 20 independent runs 
(each with a different seed for the pseudo random number generator) were made 
for each Eq and L. For each seed, at least 4 x 10 6 trial changes were made, in 
order to calculate the autocorrelation functions and their respective autocorre- 
lation times. The values we quote are the average of the values obtained for 
each seed of the pseudo random number generator and the uncertainty in r is 
the standard deviation of these 20 values. 



4. Results and Discussion 

We have simulated the case Eq = — 1 (Metropolis algorithm) as a test to our 
code. The result is presented in Fig. [5] for the magnetization autocorrelation 
time and, as we can see, all three cases have the same qualitative behavior. 
For q = 2 and 3 the dynamic exponent z takes approximately the same value 
[z ~ 2.16), while for q = 4 it assumes a higher value, namely z = 2.21 ± 0.02. 
These values are consistent with those in the literature fl6| . 

From now on we will not comment on q = 2, since this case has been treated 



in Ref. 11] . Our simulation for Eq = —0.75 is presented in Fig. [3j where the 
magnetization and energy autocorrelation times are depicted as functions of L. 
For q = 3 and L > 32 and for q — 4 and L > 64, the dynamic behavior (i.e., z) 
is the same as for the Metropolis algorithm. In Figs. [5f> (q — 3) and[3]i (q = 4) 
we present the behavior of the average cluster size < n > versus lattice size L: 
in both cases, for L > L = 64, < n > is constant. As already discussed, the 
autocorrelation function which is well described by only one exponential presents 
the greatest autocorrelation time and for, Eq = —0.75, this happens for the 
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magnetization's autocorrelation function. The dynamic exponent is calculated 
for L > L and we obtain z = 2. 16 ±0.05 and z = 2. 18 ±0.09 for q = 3 and q = 4, 
respectively. These results agree, within error bars, with the values quoted in 
the literature [l6| and with our values obtained for Eq = —1. 

For E = —0.25, our result is depicted in Fig. 2] for q — 3 ( the result q = 4 
is qualitatively the same); the behavior follows the same overall trend observed 
for Eq = —0.75, with a different value of L. For L < 128 the autocorrelation 
time for the energy is greater than for the magnetization and the average cluster 
size < n > increases with L. For L > L = 128, the picture changes and the 
autocorrelation time for the magnetization is the greater one and < n > is 
constant. The value of z is consistent with the one for the Metropolis algorithm 
(again, z is calculated for L > L). 

For other values of — 1 < Eq < we observe the same qualitative behavior. 
There is always a L, such that, for L < L, < n > increases with L and, for 
L > L, < n > is constant and the magnetization's autocorrelation time is the 
one well described by a single exponential. The exponent z, always calculated 
for L > L, is the same as for the Metropolis algorithm. This can be linked to 
the constancy of < n > in this interval: since L increases and < n > remains 
the same, the fraction < n > /L 2 decreases and the behavior is the same as a 
single-spin algorithm. 

We simulated the Wolff algorithm to perform another check of our algorithm 
and to compare with our results for E ^ 0. In Fig. [5] we present the autocor- 
relation time and < n > versus L for q = 3 and q = 4. To calculate z we used 
a different approach here 17, HI]. We perform a power- law fitting using three 
consecutive lattice size (e.g L = 512, 1024, and 2048) and call L m i n the smallest 
size. We then plot z versus l/£ TO in> as seen in Fig. |6] [l7, 18]. The power-law 
fitting for q — 3 is a good fit for all lattice sizes and we obtain z = 0.55 ± 0.02, 
which agrees with a previous estimate For q = 4, only for L > 128 we 

obtain a good fit, with z = 1.00 ± 0.02. This value is somewhat above a previ- 
ous calculation [20j for the Swendsen-Wang algorithm. If we restrict our data 
to L < 256, as in the previous reference, we obtain z = 0.94 ± 0.01, which just 
overlaps with the result of Rcf. 20] (namely, z = 0.89 ± 0.05). We plot the 
data for q = 2 for comparison with our results for q = 3 and 4: it seems clear 
that, for that value of q, the asymptotic regime has not yet been reached, for 
the lattice sizes we simulated (in fact, it is not clear if there is an asymptotic 
regime) [lH . We would like to note that a good indication of the result quality 
is the relation < n >oc L 1 I V . As we can see in Figs. 0? and[5]i, our estimates 
of 7/V from a log-log plot of < n > vs. L (namely, 7/1/ = 1.734 ± 0.001 and 
1.7498 ± 0.0009 for q — 3 and 4, respectively) are, within error bars, the same 
of the (conjectured) exact results (namely, j/u — 26/15 — 1.73333.... for q = 3 
[H and 7/1/ = 7/4 = 1.75 for q = 4 [HEl). 

In order to compare our results for different values of E Q , we have plotted z 
vs. Eq for both q = 3 and q = 4 (see Fig. [7]): we see that z is approximately 
constant for any Eq ^ and strongly decreases for Eq = 0. In fact, the subtle 
decrease of z as Eq approaches may be a crossover effect: the closer we get to 
the Wolff algorithm, the greater the value of L and one has to go to larger and 



7 



larger lattices to obtain the correct dynamic behavior. 

Therefore, we have shown that the Wolff algorithm is still the most efficient 
procedure, when compared to the generalizations for Eq < 0. But we still have 
to check the dynamic behavior for Eq > 0. We know that for Eq — > oo the 
algorithm is not ergodic. So we expect that, if a better value for Eq exists, 
when compared to Eq = 0, it is not much greater than this last value. To 
address this question we simulate the cases Eq = 0.05 and 0.1. In Fig. [8] we 
show the results for q — 3 (a) and q = 4 (b), both calculated from the energy 
autocorrelation function. As we can see, for q — 3 the autocorrelation time r 
for Eq = 0.05 is much greater than for Eq — and grows faster than for the 
latter. For q — 4, r is slightly greater for Eq = 0.05 than for Eq = 0, although 
our result is consistent with the same value of z for both cases. However, one 
has to consider that the implementation of the algorithm for Eq = 0.05 is more 
complex than for Eq = 0. 

5. Conclusion 

In this work we studied the Niedermayer algorithm applied to the two- 
dimensional Potts model with 2, 3, and 4 states. Our goal was to determine 
which value of Eq leads to the optimal algorithm, i.e., to the smallest value of 
the dynamic exponent z. We observe that for — 1 < Eq < there is a lattice size 
L, such that, for L > L, the average size of updates clusters is constant and the 
dynamic behavior of the algorithm is the same as for Metropolis'. The value of 
L increases with Eq and diverges for Eq = (Wolff algorithm). When we look 
to the auto-correlation function, we notice that for L < L the auto-correlation 
time of the energy is greater than for the magnetization and the opposite hap- 
pens for L > L. As we show in Fig. [TJ the quantity with greater autocorrelation 
time have the auto-correlation function well described by a single exponential. 
For Eq = we regain the Wolff algorithm. 

In Table 1 we summarize our findings, which show that the Wolff algorithm, 
Eq = 0, is more efficient than its generalization for Eq < 0. There is still 
the possibility that some value of Eq > may present a lower value of z, 
when compared to Eq = 0. We show that, if this value exists, it is lower than 
Eq < 0.05 and the complexity of the algorithm is greater than the improvement 
in the dynamic behavior. 



q 


-1 < Eq < 


Wolff {Eq = 0) 


2 


2.16(1) 


undefined 


3 


2.162(7) 


0.55(2) 


4 


2.21(2) 


1.00(2) 



Table 1: Values for the dynamic exponent z for the three models studied here and for 
-1 < E < 0. 
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Figure 1: Magnetization and energy autocorrelation functions versus time (in MCS) for the 
Nicdcrmaycr algorithm with Eo = —0.25 and q = 3. The main window represents the behavior 
for linear size L = 16, while the inset applies to L = 128. 



Figure 2: Log-log graph of magnetization autocorrelation time (in MCS) versus linear size L 
for the Metropolis algorithm with q = 2, 3 and 4. The quoted value for z is obtained from the 
slope of a fitted straight line for the magnetization autocorrelation time. 



Figure 3: Log-log graphs of magnetization (circle) and energy (square) autocorrelation times 
(in MCS) versus linear size L for the Niedermayer algorithm with Eq = —0.75 for a) q = 3 
and c) q = 4. Log-log graph of average cluster size < n > versus L for b) q = 3 and d) q = 4. 



Figure 4: a) Log-log graphs of magnetization (circle) and energy (square) autocorrelation 
time (in MCS) versus linear size L for the Niedermayer algorithm with Eo = —0.25 for q = 3 
. b) Log-log graph of average cluster size < n > versus L for q = 3. 



Figure 5: Log-log graphs of magnetization (circle) and energy (square) autocorrelation time 
(in MCS) versus linear size L for the Niedermayer algorithm with Eo = 0.0 (Wolff algorithm) 
for a) q = 3 and c) q = 4. Note that the values we quote for the dynamic exponent z arc those 
obtained from the greatest three values of L (see text). Log- log graph of average cluster size 
< n > versus L for b) q = 3 and d) q = 4. 



Figure 6: Semi-log graphs of z versus 1/L m ; n (see the text) for q = 2, 3 and 4. 



Figure 7: Dependence of the dynamic exponent z on Eq, for — 1 < Eo < and q = 3 (circles) 
and q = 4 (squares). 



Figure 8: Log-log graphs of energy autocorrelation time (in MCS) versus linear size L for 
the Niedermayer algorithm with Eo = (circle) and Eo = 0.05 (square) for a) q = 3 and b) 
g = 4. 
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